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ABSTEACT 


THERMAL LECOTOSITIOH OE H-BUTANE AT HIGH TEOTERATUEES 

Thermal decomposition of 9.8% n-hutane + Argon vas 
investigated in the temperature range 1060-3000K using shock 
tuhe - laser schlieren system. The density gradient profiles 
at different temperatures were got using the laser schlieren 
signals. An attempt was made to compare the experimental 
results with numerical results after assuming a detailed 
kinetic mechanism. 



CHAPTER I 


INTEODHCTIOH 

Existing studies in the literature on the thermal de- 
comnosition of n-Butane can be divided into two groups according 
to tempera t-ure . low temperature experiments (T <10001) are 
carried out in conventional reactors using pure n-Butane. 

High temperature experiments (T>1000K) are carried out in 
sliock tubes'^”'^ with dilute n-Butane (mixture of n-Butane in an 
inert gas ) . 

Thermal, decomposition of hydrocarbons is carried out in 
shock tubes because of good 'reproducibility and the conditions 
can be varied over a wide range. Many investigations have been 
done on the thermal decompGsition of hydrocarbons of low mole- 
cul.ar weight (C-j and Cg ) . The investigations on the thermal 
decomposition of n-Butane is incomplete. There are differenres^ 
among the results at lower temperatures and the ones at higher 
temperatures. There is a marked difference in the data for the 
order of overall react! on, and there is a disagreement on the 
product distribution also. Some authors ’ found the order of 
overall reaction to be 1 while Sagert and laidler and Blafehmore, 
Barker, bud-'Ccr-an'oran'^ postulate the order to be 3/2. 
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At higher temperatiir es there exist three studies. 

First one is hy ¥ittig^, using a single pulse shock tube, the 
second one is by Izod, Kistiakowsky and Matsuda , and the third 
one is by Fritz and Gronig'^ using a shock tube coupled to a 
tirae-of-f light mass spectrometer. 

In the present "work, n-C^H^Q + At mixtures containing 

9.&A n-Butane ■were heated by incident shock waves over a 

temperature range 1060 - 3000K using shock tube laser-schii=iren 

R 

technique of Kiefer and Lutz . 



CHAPm II 


EXPEEIMENTAL SETUP ASL LASER SCHIIEREN SISTEIi 
Experimental apparatus : 

9 

The experimental setup is completely described elsewhere . 
A diagramatic sketch of the experimental setup is shown in 
figure 1 . 

Th e shock tube ; 

The driver section consists of a 15*3 cm I.E. 1.68 m 
long stainless steel cylindrical tube. The low’ pressure section 
Q.e. the driven section or the test section) consists of a 
9.85 cm I.E. 5.72 m long honed stainless steel cylindrical tube. 
The driver and driven sections are joined by an intermediate 
tapered stainless steel section which fits smoothly with the 
driver section on one side and with the driven section on the 
other side. This tapered section is permanently attached to the 
driver section. The driver and the driven sections can be 
locked together or unlocked by means of a coupling device. The 
entire assembly is mounted on wheeled iron stands w;hich in turn 
are mounted on iron rails.. 

The driven section can be evacuated to a pressure of 
2 X 10”^ torr using a 4" silicone oil pump backed by a 
roughing pump. The driver section is evacuated using a separate 



4 


mechanical p-ump. The comhined leak and degassing rate in the 

driven section has heen measured intermittently and found to 

-S 

he less than 5 x 10 torr. The pressure in the driver section 
is measured using a bourdon pressure gauge. The pressure in 
the driven section is measured using a silicone oil manometer, 
-hock v-’eves are generated by bursting mylar sheets of thickness 
0.005" by a pneumatically operated mechanical plunger which 
runs along the axis of the driver section. Mixtures of H 2 ,N 2 , 

He and 00^ were used in the high pressure section to obtain 
the required variations in the shock velocities. 

The shock speed was measured by means of gold resistance 
thermal gauges mounted flush with the walls of the shock tube. 
Since the shock wave attenuation was found to be negligible 
only two resistance gauiges were used to obtain the shock 
velocity. The thin film gauges were prepared by vacuum depo- 
sition of gold on to glass tubes closed at one end. The dia 
glass tubes were sealed at one end after inserting platinum 
wires. ¥hpn this end was polished the cross section of plati- 
num wirps appear as small circles. Gold was then evaporated 
under vacuum on to this end to form a thin film of 1mm width 
making electrical contact between the two platinum leads. 

The resistance of the thin film plate was continuously 
monitered during vacuum deposition and deposition was stopped 
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■when the resistance reached 50-/V. The thin films'were mounted 
on the shock tube with the longer dimension perpendicular to 
the axis of the shock tube. The first gauge was placed at 
4.25 m from the diaphragm and the second gauge was placed at 
0.785 m from the first gauge. The windo'ws for passing the laser 
b'='am were placed 5 cm after the second gauge. The resistance of 
the thin film gauges changes due to the passage of shock across 
them and it was converted to voltage signals which were suit- 
ably amplified and used to start and stop a Beckmann model 7370 
digital counter. The time recorded by the digital counter was 
used to measure the shock speed. The rise time of the thin 
film gauges was found to be about 1 psec. 

G ases : 

High purity Ifetheson n-Butane (99.995^°) aud ultra high 
purity Argon (99-999%) from Indian Oxygen Ltd., were used. 

Mixtures of n-Butane and Argon were prepared in a 
thoroughly evacuated stainless steel mixing chamber. The gasses 
arp allowed to stand in the mixing chamber for 24 hours before 
use. 


Laser Schlieren System : 

The observation station i.e. the point where the laser be; 
'©5‘®8 through the shock tube was located 5 cm downstrecam of the 
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spcond film gaup:e and 3.5 cm Isefore the end plate. The reflec- 
ted shock vvave did not interfere v,ith the measurements because 
they (relaxation process) "were usually completed within lOpsec. 
behind the incident shock whereas the reflected shock would 
arrive at the observation window? only after 50-100 psec. 

Class windows of 1" dia in plexiglass adaptors were 
mounted flush with the wall of the shock tube. The windows were 
cleaned regularly with the lens paper to ensure gaussian 

O 

distribution in the transmitted laser beam. The 6528 A laser 

beam from a 0.5 mw Spectra Thysics ilodel 156 He-Ne laser passed 

through the shock tube through these windows normal to the axis 

of the shock tube. Pepeated runs in Argon, with the laser 

removed and realligned did not show any difference in the 

resolution, showing that the perpendicularity was good. The laser 

beam coming out of the shock tube was reflected by an aluminum 

mirror to fall on a knife edge. The distance from the shock 

tube centre to reflecting mirror was 1 .35 m and the distance 

, 1 . 

from the reflecting mirror to knife edge was 5.49 m. The part of 
thelight that was not cut off by the knife edge was collected 
by a 7 cm dia lens and focussed to a H.P. type 5082-4203 pin- 
photodiode. The pin-photodiode along with the associated 
amplifier was mounted on a device which could be adjusted such 
that the photodiode receives all the light collected by th»-'lens. 
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The knife edge covld he moved smoothly to cut the laser team. 

The pin-photodiode assembly, the glass lens and the knife edge 
■were mounted on a smooth rigid bench. The laser and reflecting 

'''^f're mounted on separate stands. The reflecting mirror 
was rotated at a slow uniform, speed so that laser beam sweeps 
the knife edge at a unifoi^m rate. The voltage output from the 
pin-photodiode as the beam sweeps across it is shown in figure 2 
It is seen that the curve is linear in the centre and linearity 
extends to + 209-0 from the centre. 

The photodiode output, suitably amplified, was fed 
throueh a Tek type 1A5 plug-in amplifier to a Tek type 549 
storag'"'’ oscilloscope. The oscilloscope was operated in a single 
sweep mode "which was triggered by the signal from the thin film 
gauge just ahead of the observation station. 

The overall response time of the systom was measured 
by firing shocks into Argon. It was found that schiieren record 
decays with a time constant of 0.12psec. which may be taken as 
the time constant for the detection system. So schiieren 
records which showed time constants less than 0.35 p 
not used in analysis. 


sec. were 
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Pr ocpcj urr^ : 

Thp driver end driven sections i^ere locked together by 

means of a mylar diaphragm in between them. The driven section 

was evacuated to about 5 microns with the help of the two 

roughing pumps. The pressure was further brought down to 
—5 

about 5x10 torr using the diffusion pump. Runs taken with 
evacuation times varying from 2-10 hrs . did not show any 
difference and hence the evacuation was usually can ied out 
for 2 hours. The driver section was evacuated using a roughing 
pump and was then filled with the gases of required compo- 
sition and pressure. The knife edge was moved away from the 
lens arid a chopper was placed in the path of the laser beam. 
This produced an a.c. signal from the pin-photodiode with its 
peak corresponding to full beam falling on the photodiode. 

Then the knife edge was brought back slowly cutting the light 
falling on the photodiode such that the peak voltage was equal 
half thefull beam voltage. The beam was thus centered and the 
chopper removed from the path of the beam. The tost gas 
(i.e. mixture of n-Butane and Argon) was now filled into the 
driven section to the desired pressure and the shock was 
fired by releasing the mechanical plunger. The counter time 
was noted for shock speed calculations and the schneren 
record stored in the oscilloscope was photographed. 



CHAPTER III 


experimental results* 

The shock speed ^^as calculated "by using the distance 
"between the two film gauges and the counter time. The equi- 
librium values of the temperature, pressure, density etc "were 
calculated for each shock speed by solving the mass, momentum 
and energy conservation relations by an iterative procedure 
on a computer. The mass, momentum and energy conservation 
relations are given in Gaydon*'^. A listing of the computer 
urogram used is given in Appendix 1 . The specific heat data 
is taken from NASA SP 3001 . The equilibrium values for five 
selected shock speeds are given in Table 1 . 

Lensity gradient profiles behind incident shock waves 
vv'tn e obtained in 9.8% n-Butanet- Argon mixture for five diffe- 
rent temperatures in the temperature range 1060 - 3000K. The 
oscilloscope traces of the laser schlieren signals for the 
five different temperatures are shown in figures 2-4. Table 2 
summarizes the range of experimental conditions for 9.8% 
butane + Argon mixture. 

Lensity gradient profiles can be obtained numerically 
by solving the mass, momentum and energy relations, after 
assuming the detailed kinetic mechanism. The detailed kinetic 
mechanism is discussed in the next chapter. 



GHiPTER IV 


KINETIC MECHANISM 

The pyrolysis of hydrocarbons are recognised as free 

radical chain reactions. Rice and Herzfeld"'^ first proposed the 

free radical mechanism in discussing these reactions. Sagert 

and Laidler^, Wang, Rinker, and Cor-coran^ and purnell and 
1 1 

Quinn have proposed reaction schemes for the thermal de- 
composition of n-Butane, based on the Rice and Herzfeld type 
mechanism. 

Initiation of chains in the pyrolysis of n-Butane is 
by the decomposition of n-Butane into tv/o free radicals, 
decomposition of n-Butane primarily takes place in chain 
propagation steps. The chain termination step of the mecha- 
nism is the blraolecular recombination of two free radicals. 

Sagert and Laidler^ have proposed an eight elementary 
step mechanism for the temperature range 793 - 863^. Wang 
etal^ have presented a 19 reaction mechanism for the tempe- 
rature range 733 - 833K. They have considered the formation 
of propane, Butylene and the propagation of propyl radical. 
Blakemore etal"^ have investigated a 13 reaction mechanism in 
the temperature range 803 - 873K. They have encountered the 
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formation of Butylene and fhe propagation of propyl radical, 
but they have not considered the formation of propane. Yexj 
few investigations have been done at high temperatures. Pritz ' 
and Gi'onig^ have investigated the thermal decomposition of 
n-Butane in the temperature range 1 1 70 - 1570K. They have 
coupled a shock tube to a time -of -flight mass spectrometer. 
They have developed a 21 reaction mechanism involving propane. 
Butene and propyl radical. V'ittig^ has investigated the 
pyrolysis of n-Butane at 1400K using a single pulse shock 

tube. He has used a mechanism similar to that of Sagert and 

2 '5 11 

Laidler , Wang, Rinkbrr and Con'r^coran and Prunell and Quinn 

but modified for the conditions of high temperatures. 

A free radical mechanism involving 14 elementary 
reactions is given below. 




CH^ + C^Hg 


(3) 
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Chjijn 




k, 


° 2®5 ^ ‘^ 2^4 


CHj + Q L_j O^Hg + OH^ 


C.H 


2--5 


-» CoH, + H 


2 4 


H + 1 » o^Hg + Hg 


2CH. 


k 


8 


C^Hg 


Cr^Hp. 4“ H 

2 5 


k. 




C„H^ + H 
2 5 


k 


10 


-» H 2 + CgH^ 


CgHg + H 


k 


11 


•-*. GJic+ 

2 ? 2 


CH^ + H 


k 


12 


CH_ + H„ 
? 2 


C„H^+ H 
2 0 


•13 


-fr 2CH. 


Termination: 


202Hg 


11-* O^H^q(OE) OgHg + 


(4) 

(5) 

( 6 ) 

(7) 

( 8 ) 

(9) 

( 10 ) 

( 11 ) 

(12) 

(13) 


(14) 



15 


In the mechanism given above the main products are 
ethylene, propylene, ethane and methane. Besides these main 
products there are reaction products such as hydrogen, traces 
of propane, traces of i-butane and components ■with higher 
carbon content. J\cetylene and carbon were formed only at 
higher temperatures and longer reaction times. So they were 
not considered in the reaction mechanism. 

The frequency factors, activation energies and the 
rate constants at 127'5K for the above mechanism are given in 
Ta ble 3 . 


In the mechanism given above there are 10 species. 

Pi mass balance is made for each component and the corres- 
ponding differential equations are given below. These 
equations are differential equations of concentration with 
I'cspect to time. 


- k2 [ CgEp - jpii, ] 

-1^7 [ Viol [S 1+ ^14 [Vsl ^ ^■'5^ 

^ 4 ^ = IC 2 [C^H^q - k^ [C^Hg] - k^ [ ] 

+ kg [O^H^O ][CH3] +k7 [C^H^o][H] (16) 
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dt 

d_[C2Hg ] 

dt 


d [CjH^] 
dt 


dt' 


i_ ^^^’4 ^ 

dt 

a [ CH^ ] 
at 


a_[ Hg] 

at 


= k3 [o^Hg] (17) 

= kg [C^H^q] [CgHg] + kg [CHgll 2 +kg[C2H3][H] 

- k^^r.CgHg] [H ]+ [CgHg]^ (18) 

= - kg [C^H^g] [CgHg ] + k^ [ C^Hg ] 

- -^9 m -'"10 P2%]F^ 

+ k,^[C2Hg][H] - k„ [CgHglEH]- 2 k^^ [CgHg] ^ 

= [C^Hg ]+ kg [OgHg] + k^gl OgHgH H] , 


tk^^ECgHgl^ (20) 

= kg [o^H^glECHg] -k^gEcH^lEH] (21) 

= k; [04Hg ]- kgE C^H^glECHg] -kg [CHg ] ^ 

+ k^j[CH^] [H] + 2k^3[C2Hgl p] (22) 

= k^ [K] + k,Q[ CgHgl [H] +k^^ [OgHglEH] 
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- [C^HglLH] - k^2 PH^] [H] - k^5 L C2H5] [H] (24 

The 10 first orfler differential equations (15) - (24) 
are solved sim-ultaneously using numerical methods. The initial 
conditions for solving the above differential equations are 
initial concentration of the respective species. In this 
case all the initial concentrations except that of n-Butane 
are 0 . After solving the above equations the concentration 
profiles of the different species can be dravin as sho'wn in 
figures 5 - 7 . 


iJ H] 

dt 



CBAPTER V 


NUMERICAL RESULTS AND DISCUSSION 

A 19 reaction mechanism involving h-utylene, propane 
and propyl radical was also investigated. The reactions that 
were added to the earlier mentioned 14 reactions are as 
follows 





CH, + OjH^ 

C4.H9 

^ 

O4H8 “ 


^ 

CjHg + H 

C^H? 

^ 

+ CHj 

CH^ + C2H5 - 

^ 

SB 


The formation of butylene, propane and propyl radical was 
very small ( of the order 10”^*^ ) and it continued with 
longer times. So for the final computations the above five 
reactions have not been considered and the 14 reaction mecha- 
nism which has been discussed earlier is used. 

The simultaneous differential equations are solved 
by the Hamming's predictor - Corrector method. The algorithms 
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1 2 

and computer program are discribed elsewhere . A listing 
of the program used is given in Appendix 2. The initial 
points for the Hamming's method are obtained by applying the 
Eunge - Kutta method over the first three steps. Then the 
equations are solved over each step by applying the predictor 
and corrector algorithms alternately. 

Figure 5 shows the concentration profiles for the 
decomposition of n-Butane and the formation of butyl, ethyl, 
methyl and hydrogen radical at 1273K. Figure 6 show/s the 
concentration profiles for the formation of propylene, 
ethane, ethylene, methane and hydrogen at 1273K. Figure 7 
shows the concentration for all the species at 1391K. Con- 
centration computations were made till 2742K. At higher 
temperatures the decomposition of n-Eutane is so fast that 
the integration step size for the time had to be given in 
the range 0.001 4 sec. 

To make a comparision with the experimental results 
the differential equations involving the density and the 
temperature have to be- coupled with the equations (l5) - (24) 
The initial conditions for the density and the temperature 
are the equilibrium values (i.e. the frozen shock conditions) 
But unfortunately due to lack of time the last part could not 
be completed and is left for the future. 
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NOMENCIATURE 
Initial density (g/cc) 

D 2 Equililrinra density or density at frozen shock 

conditions (g/cc) 

^21 ^2/^1 

kj^ Rate constant (sec"'^ or cm^ mole sec ”* ) , i=1 ,14 

Initial pressirre (torr or ram Hg) 

Pg Equilihrinm pressure or pressure at frozen shock 

conditions ( Of ) 

P21 P2/P1 

Ipjl Temperature at frozen shock conditions or equilibrium 

temperature^® K j 
u^. Shock speed (mm/ P- sec) 

0 

rc^i'I 1 Concentration of C^H (moles cm 

n n m 

n = 0,4 and m = 1,10 
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TABLE 1 

EQllIIIBEIUM VALUES EOR PIVE SELECTED ’ SHOCK SPEEDS 


mm/ u sec . 

s 


p.^mm Hg 


®21 

1 .027 

1273.23 

2.88/ 

■ V 

12.8521 ' 

3.0095 

1 .135 

134-1. 33 

3.92 

16.5886 

3.5548 

1 .596 

1972.48 

6.08 

36.0070 

5 .4427 

1 .877 

2395.23 

7.80 

50.9941 

6.3474 

2 .091 

2742.14 

1 .68 

64.0433 

6.9634 


TABLE 2 ■ ■ 

RANGES 01’ ' CONDITIONS EOR 9.8°/ n-C^H^^ + Ar MIXTURE, 
ram Hg Ug mm/ Usee. T^ K 


1 .52 - 7.80 


0.849 - 2.302 


1060 - 309® 


TABLE 3 
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ASSUfCED MECHANISM AND RATE CONSTANT EXPRESSIONS USED IN 

COMPUTED SIMULATIONS 


Reaction Frequencj?’ Activation Pate constant Refe- 

factor* energy at 1275'^* rence 

Kcal mole 


1 

1 .0 

X 

10^7 

80.0 

1 .85 

X 

10 ^ 

2 

2 

7.7 

X 

lo ”* ^ 

10.4 

1 .26 

X 

10^0 

2 

3 

6.5 

X 

10 ^' 

24.0 

4.93 

X 

10 ^ 

2 

4 

1 .6 

X 

10 *'° 

22.0 

2.68 

X 

10 ® 

2 

5 

2.7 

X 

lo '*’^ 

9.0 

7.70 

X 

10 ^' 

2 

6 

3.0 

X 

10 ^^ 

39.5 

4.97 

X 

10 ® 

2 

7 

1.0 

X 

10 *''^ 

7.9 

4.40 

X 

io '® 

2 

8 

1 .0 

X 

10 ^ ^ 

0.0 

1 .00 

X 

lo ’’’' 

3 

9 

1 .0 

X 

10^5 

0.0 

1 .00 

X 

10 ^® 

!3 

10 

1 .7 

X 

lo ’' ^ 

0.0 

1 .70 

X 

lo ''^ 

13 

11 

1.3 

X 

10 ^^ 

9.4 

3.25 

X 

lo ’'^ 

13 

12 

7.2 

X 

lO ^'^ 

15.1 

1 .89 

X 

10''2 

13 

13 

3.7 

X 

10^5 

0.0 

3.71 

X 

10^5 

13 

14 

1 .6 

X 

10^5 

0.0 

1 .60 

X 

10^2 

2 


— 1 5 —1 -1 

* In the units sec” or cm"'^ mole sec 
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Echerr'atic diagranr> of shock tube. 









2 -Loser Schlieren sio.nal and the corresponding 
calibration .urve(i.c photodiode output as 

laser beam sweeps across the knife edge) 
Pp1-68torr , us r 2-2492 mrr,/MS, T = 274 2°K 
Laiibrntion trace - Horizontal scale r 0- 5us /div 

Vertical scale 

Curve a =0-5 Volts/div. 
Curve b =0- 2 Volts/ di v, 





■ Loser Schlieren signal in 9'8°/o C 4 hlio-"Ar 
P]=4-32torr, us = 10304 nnnn//JS,T2=1276°K 


+W+4-|,H4f-4-H+4U 




Loser Schlicren signal in 9-8 /o C 4 Hio Ar 
pi =3-92 to rr , us = 1-1353 mm//JS,T2 =1391°K , 







++-1- 1- +-1-H-H "tft-H ■H+f-ht+ii+r 


"**1 r~2yus 

Fig. 4 a-Las 0 r Schlieren signal in 9-8 % C 4 Hio + Ar. 
Pir 6 -O 8 torr, us=: 1-5965mm//js , T= 1972 °K. 


-Laser Schlieren signal in 9-8 loCi^H]Q*Ar 
Pl=7-80torr , us =1-8779nnm l^s ,T= 2395°K, 










Concentration , moles cm 
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Fig. 7 -Concentration profile for p2/p'i = 16-5, T= 1391®K and Pi = 1 torr 
C4Hio:Ar = 9-8:90-2. 

10^CC4Hio],10’’ CC4H9],10'^CC3H6],10’[C2H63,10”[C6H5] 
10’^[C2H4J, 10^ [CH4], 10^[CH3],10'7 [H],10® [H] 




J^PPiCMDTX 1 


" r'-ifc; PRD'^K^’^ ro nHpiJfE 

iT'ir-trRft r^iR^; ^^H•) DE'iSITf 


T!}E E3UIl.tRRXUM VALUES OF 


’ZM2X,'CP?1',12X,'P21M4X,' 
SHOCK calculation WITH BUTANE 


A LISTTv: 

PRESSURE 
ILINsSO 
TZ = 0 

0 = l2PbnO.U „ 

OPEN (UNI r=1 jFTLEa'SUnCK.CDR'I 
TYPE^On 
TyPE200 
DO 502 lal,25 
READC IrSnon.Ul ,n21A 
FORMATCiFlO.S) ^ ^ 

rORMAThx,rl5.5,2X{6CF14,8,2X)//5 , 

FOR«ATC2X;//l5x,'Ul%12X,'T2'r10X,‘ 
lD2t',13X.'y'///) 

rORMATClHO, 10X,70HRUTANE+AR FROZEN 

1 CONCENTRATION a 0,09800///) 

lla8,317E + 07 • 

Rial, 987 

Tl*298.15 

Wl»41,724 

X*0, 09800 

CP1»5,2S20 

Alla+3,276880Be0O 

A21»+7,080814E-01 

A31**5,833ll8E-06 

A41»-1.5832666E-08 

A5l«0,6 

ASlaOjO 

811b+2,P43T494F00 

R21S+8.31RS78E-02 

R3la-3.R^ 31026E-P5 

B41s+7,i:»11 3E-0 9 

351*0.0 

361=0,0 

Wll=5fl,0 

A12=+3,l364326F0O 

A22=+l , 3(>HOhORF-03 

A32=+2,l837i64E-06 

A42=-3,y518236E-12 

a52=0,0 

A62s0.n 

B12=+3,9USiSlEt>0 
322=+7.a‘>66J 07F-04 
33?=-3.2'’3Sl 77P-.07 
342=-»'6.3Sl 32Q 3P-1 I 
352=0.^ 

362=0,0 

SVI 2 = 48,1)66 

A13=+3,O2l2bO1F0O 

A23=-2.1737249E-03 

A33=+3.7S4220 1F-06 

A43=-2.9O47200F-09 

A 5 3 = 0 , 0 

A63«0,0 

B13=+2,5372567E0O 

R23»«l,9422l90E-05 

Sl3*-8,«017§21E-O9 

B4Sa+5,9543C>2lE-12 

BS3«0,0 

B63»o!o 

W13»16,00 

A14e4.9675 

A24*0.0 

A34=0, 0 

A44a0,0 

AS4sO, 0 

A64=*i,4S10601E03 

Bl4a4,9765 

»i4*0,0 

B34fc6,0 

844*0.0 

ll4*2i!48l0601E03 

Hl»CCAii4n+A21*Tl4»2/2,.fA314(T14*3/3,'l'A414fl4*4/4,)*X4(Al44'Tl+'A6 

i4/Tl)>T(J .-X) ) 

xi*xni,-Y>/ci.Ty^x) 

X3«T»X/Cl,H-Xf:X) 

X4»£l.*X)/(l.+y*X) 


10 


12 

13 

14 

16 

17 


la 

102 

9 

02 


41aftU»KlfAi2*X24-Al3»X3 + Al4>|tX4 
A2sA21<'Xl+A22*K2 + A23»X3 + A24>l'X4 
A3sA31»XltA32«'X2 + A33>l'X3 + A34*X4 
A4sA41»XI+A42*X2+A43*X3+A44*X4 
A5=A51#XlfA52tX2+A53*X3+A54*X4 
A6aA61*Xlf A62*X2+A63»X3+A64^X4 
R1«B11*X1+012>»«X2 + SJ3*X3 + 014>I'X4 
B2a821«‘XlfB22*X2 + 823*X 3 + 0244X4 
83*B3 I #XJ +0324X2 + 83 3*X3 + B34 4X4 
84=04 1*X1+0424X2+B43*X3+B444X4 
B5sB51*X1+R524X2+B534X3+B544X4 
B6=B6l*XJ <-0624X2 + 8634X3 + 0644X4 
W2*l4114X3 +W124X2 + W1 34X3 + W144X4 
P2t»021A 
P12ai,/D21 

P2lal,-+CWl*Ul**24(U-Dl2)/CR4rJ)) 

DH2A»€tll442}*Cl ,-(012442) )/8,36E+'07 

T21»Ph4Dl24W2XWl 

T2»T2l4ri 

Z«T2*4C-l,/3,) 

]EFCT2**790 *) 12 J2 13 

H2r»Al+A2lT2/§,+^34T2442/3,+A44T2443/4.+A54r24*4/5, 
CPT2=A1 +A2*T2+A34T24*2+A44T2*43+A54r2444 
no TO 14 

82Ts0l+824T2/2.+R34T2442/3,+B4*I24*3/4,+B54r2444/5, 
CPT2=B1. + 0 24C2/?,+834T2*42+844r7.44 3 + B54r2444 
H2=CfH2I*r2-'-U )/W2)+CY*X4D/CCl,+y4X)4M2}) 
CP2t.=CPr2/CPl 
EPS*0,0001*t)^»2'' 

TF(ABSfH2«nH2A)-EPR)17,17,16 

021=021 +4,^>riOf'i 
nu TO 10 

T jj 2 e j 7 4. 1 

TF(iX-TZ/U,r'v4Tt,iN)i9 18 19 

TYPE! 02 

FORMATfl H J ,///) 

Typ(.;t 00,!I1 , r2,7,,CP2l,P2l,D2l,Y 

rn^TTNOF 

STOP 

Kf-f) 


+A6/T2 

+B6/T2 
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C’ A PKUUKA'I I'J SULVe blMULl’AfJKQUS Dit’FEKeNTAl. EaUATlDNS 


C 

C 

c 

c 

c 

c 




»i'Hi.s pro-iram solves a system of 
»ntial equations using tne rfaining 
»metnoci, me program reads a star 
♦varlaoie x,tne integration step 
^integration XMAX.and n initial c 
»are macrices containing solution 
»aescrlb»a in tne text.TECJ) is t 
»tor tne corrector eguatlon.coUNi' 
»rUNCXiuM KUNGE l8 called to inte 
♦steps and yields tne starting va 
♦method; mere after the equations 
♦eiy tne predictor and corrector 
♦maminG, me equation eeing solved 
♦,ytoj»u,y,dycyJ/dx»i,u,to solve 
♦equations ,aii -symoois refering 
♦oy statements detlning tne cieriv 
♦ the solutions are pRiNi’ed at tne 
♦iNt integration steps there afte 


N first order ordinary difteref 
•6 predictor - corrector * 

ting value for tne independent* 
sl«,e H,tne upper limit of * 
ondltions XRtt)..,YRCW3,'y & E ♦ 
and derivative values and are* 
ne truncation error estimate f 
is tne step counter, rne * 

grate across tne first three * 
lue needed for tne Hamlng's * 
are integrated using alternat* 
portions Of tne rUNCfioN ■ * 
as an example are .ozy/DXZj-y f 
another system of dlf f erentlalf 
to F and FR must be replaced ♦ 
atives for the new gystem, ♦ 
initial X, and after every * 
r • ^ 

y; y ;K m f y >te 3|; f !|( f f m m f )|c !(! f it; f f f 3|C !|( jfc 4^ 


1 


liMTEGtH C;jUNT,RU'VGE,HAMiNG 
HEAP MrK/',K3,K4,K5,Kbm7,Ktt,K9,K10,Kl 
Oli^KNaJ, JiW TEC20),lfRU05,FR(2U) ,XC4#20) 
ICIU) 

UPE.nuNl J' = l ,DE VICE®" DSK' e FILE® 'ASS, CPR 

HEAPCl,'^) MAX, INI’, N, ( XR C J ) » d«l » « 3 # 

■riPE 2UU, H,XMAX,iNT,N,U,Jsl,N3 

TIHE 2U1, X, CXHCJ) f 

Kl=b,U13E'd 

R2 = .8t)otl 1 

K3«,eiy /ElU 

Kf»,ib7iE¥ 

KJ>»,4U7»eil 
Kb*,74blEH 
K7*,lSU2E14 
Ka*i,6Ei J 
xy=i,oEi j 
K10B1,UE11 
Ml»l, /E+12 
M2al,«4i2Eli 
Mi»J,UObEli 
K14®J, VlbCl J 
r«82.yb 
rr»i,yB/ 


l,Kl2,K13,Kl4f HIXMW 
,rCi,203 ,W(105 ,CPRC103 

'3 

tMW(J3 ,J=l,N3 


2 

i 


R=l/ /bO.U 
CUUNTsU 
UU 2 tJ = i, 
xE(;ij3 ■» u 
yt4,j3 « 
irCRUNGEt 
r«(i3s-xj 

llOJ+KH^t 
rKC2)=(R2 
irR(l3 3»7 
FKCi)aK31« 
rK(43«(K2 
i2+KH»yKC 
r H C E> 3 S ( 2 , 
lRCb3»»2- 
ltS3^XRtl 
rRtb3aK4» 
FKt7)»(K5 
rRt83«K3» 

FHtS3*(R7 

ItlOS^YRC 

F«{103»K6 


N 


IKt J3 

N,rR,rR,X,H3,Ne,13 go 

♦ XR(i)-CK2^yK(;i3 3^XRCE> 
lCRtb3)*^2 

fYK(133»XRt53-CK3^yRC2 

Rtl03 

VK123 

♦YR(1) )*yKC53+K9»YRClU 
53*2-K12^yRC43^YR(l03 
♦K1 )^YRCi3-CK2^YRCl)3» 
K9yyRClU3*YRC53-Kll»YR 
U) 

YRC 2)fKo^YHt5)+KR^YRC 
fYRtl33^YRC8)-K13^yRC7 
YKU3’'CKb^XRtX33*YR(8) 
53^YKC103 

♦ YRC133*YK{.10)+Xli^YRC 
/3 

1tYK(53-Cl«l7^YRC13 3»YRtl 


)2cKb^YKU3 3»yR(83-CK7^YRCl)3*l'RC 
33-K4»YRC2) + CK5fYRa3)*YRC8) + (K7* 

)^yRU3+KlO^YRC83'' 

yKCb)-«-K4^yRC2)-K6?yRC5)*»2,*CKt)*Y 
(;5J^yRtl) + Ki2*YR(4)tYRti03**Ki4*yR 

^ «. fW 4, 

53'*2+Kll^YR(5)*rRC10) 

)*yK(l03 

«Kiy^YR(«)'*2+Kl3*yR(7)^YRCi03t2,. 

b3*YRU03tK12^yRC43^YRCl0)+Kl3*yR 

0>-K9fyRClO)fYRC53*KllSYR(5)^yRCl 


CUUf^T 3 CJUM!’ f 1 
ISUH 3 4-'JJ‘.INl' 
i^U ^1 3 1 f ^ 

¥l ISUi*, J) sYKl J| 

UM*?HU *)*YHUUJ ^ ^•'^'^♦’f«l 2 ) + CK5»yR(i3 )*yRC 8 ) + 

KtlSUH, J J 3M*YKr^5 

4l » i • ^ ♦ K « » t H O ) ^ Z - M it » Y K C 4 3 ♦ Y K U U J 

l43YMtlUl*YKtt>> * • ... ... 

tf i l»U»,t»3 8«.%f JKlifJtRtifYKl &JfK8*xKl5)-;^ + Ml»XRt53 >|t3(:R(,lU) 

rtiauts, /;»tl\!»»*«U33*YKl«3-K14»YRaU)**KCn • ■ 

r l l»UH,1» J «IK/»IKU) 3♦^MtlU)+Kll*YK^5J»^KClU J + K12>''XK(43»XRCiy) + Kl 

lJi*IHUW3 *»Kl 3 J • . 

rtiaua, lU}»^tt*lfK(.!>|•l K7»YKCU)*XKUUJ-K9>t'yRClQ)^2fRt5)-Kll’t=XR(5) 
J'«‘!EKUU;-M4J*IKi4j*1tKUO)-Kj3’l«KRt n»XRUO)-Kl4*YRt53*)fRUO) ' 

IP icuuNivti'ii^twr, "K.cntirn*) nrj to 7 • 

ir tCUUM'.-ljK.i) iit*h jt, t7USUB,a3 , Jsl.ftf) 
ir ccoMwr ,!S'fi 33 tyrb; 2 y;^rX, cyu WJ j > m> 

JLrtX*«r*X'1AJl«H/2, 3 GO TU 1 
4 riCOONi.Gr.i 3 GO fO J 
M»HAMlNGtN, I ,rrX,H-iE3 

rUf 13**K4rit(lrl3’*K2*iCl,UrtUf&3-K&»iUrl)^y(lf »)'-K7*1fCi,J.)ryC 
Itf lyj+aarcptifSi)**;^) - . ^ . . 

rtii23»tR2*m.l3ryCl,S)*K3*ycn23-K4»iCU23fK5ryci,13ryCl,8)+K7 
trYti#i3*ut,iy3 . . 

rtw ^5=<3>¥a.a3 

rtl,4jsi^^^f*(i,l3 3r3f(lfS>3+K‘#*XtlflOJ*Ytl»53tKlO»yc 
ril,*>J«iJ.*!^l»Xtl,13»K2rXtIrl3»yCl,!)3'«'R4»y(l,2)-K&ryclrS3*"2,»K8»y 

Ul I l,iy3*yu #i> 3 -Ml’’'y u ,b3*y 13 ;10)+Kl2»5fcL45ryUf 103 

J*M487lS,^3*Ytl,lt33 ■ ■ . . 

rtl ,«13 3 K <»»7 t» , 23 +«¥b»Vtl fb3+p.84Xll,b3‘‘2+Ml*yCl,5)8yti,l0) 

p Ur /3»l«.*>»ttl,t33*ttlf83-Kl3'('Xtlr73^XU,i03^ 

PCI, w 3 SR <»x C J ,23-KS>'PCXll,l )=PYa,'«3 3«-MU»y U,83'‘2+Kliryci,73*ya, 

PuJJisiK/M 11^33*4 til iy 3 +Ml»Ytl,b 3 ’PYCl,lU 3 +Kl 2 »yti, 43 » 3 rCl ,103 
l+M3*Y 1 1 ; / 3*Y Cl , iy3 ...- 

rci,iy3sK;*)((i,b)-tK7»Yci,i))»yti,io3-Ky»YU.io3na.a3-Ki|»5fti| 
ia)¥¥ll,l03-Kl2*YCl,43»Ytl,4y)"'Hl3'PYCi» / ) ♦¥( 1 , 1 0) -Kl IfY 1 1 , 5 3 ♦¥ C I 
1/103' ' ... 

sauMisy.o 

DU bO J»t,N 
aauMi»ssu«i +YC 1 ,0) 

S»«IM2W0,0 
»0 bO J*l,N 

sau«2»sa!»N2 + CNNtoj»i[ci,a3 3 

MlXffWttSSUN^/.SSUNl 

TYPi; 

irtM.isa.i) GO TO n 
COUNt » CJUNf + 1 

P’UKMAi Cl>X,P'iO«4,3lOX,FXO,6fl2X,riy.4/5X,i5,15X,I2/ ,,,^ 

%%f7hHU/^lHyj5xiJlHx7l&Xr4C2HYC?I2j?H3,iyX 

.13/ tlH /tbX, 4C2HyC|i2,.lrt)»iyX)3 )_ ^ 


ij/ tin flOJVr .p,* ,, ^ 

ruKMAft'/ 4 H , Bib;/, 4i:ib;7/ClH ,iyX,4SlS,73 3 


BK13 
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t UN CTi UN fl AM I NG t N , y , r , X , H , TE ) 

♦ HAMiNG INPUKMSNl'S MAMING'S fc'KEUiC t’Urt-CORRECI'OR AGUGORiTHM TO » 

♦ SUGVE W SlMUU'i’ANEOUS FIKiST UKDEK PirfEHCNTlAG EUUA I'lUNSCORDT) ♦ 

*X IS THE » 'JDEHE’i'JE rr VAKiAttt.E ANU H IS TriE INTEGRATION- STEP * 
»&liiE.rME HJIU'lNK MS! oE CAGGEO TWICE FOR INIgGRATlON ACROSS » 
»KACH STEi^.J;! TME KIKST CAlib IT IS ASSUMED THAT IME SOOUIXQN » 
»VAtiUES AJ‘1 THE OEHIVATAVt VALiUSS K'JK THE H EJUAliUNS ARE » 

*STUKED IN I'HE FlKST N CUDUMNS Jf T.RE FIRST FOUR HONS OF THE ♦ 
*X MATRIX AND FIRST THREE HUNS UF THE F MATRIX KESPECTIVELT , » 

»THE ROUTINE COMPUTES N PREUlCrEO SJUUTIONS TPREDUJ , I NCREMRNTS* 

♦ X ST'M'APIt) PUSHES ALU VALUES IN X&F MATRICES DOWN ONE ROW, » 
MTHE PREOtCTEU SULNS IPkEDIU) ARE M 00| FlED , US ING THE TRUNCATIONN 
♦ERROR ESTIMATES TEU) FROM THE PREVIOUS STEP, AND SAVED IN THE ♦ 

♦ FIRST RUN OF THE I ' MATKI X , HA MI NG RETURNS TO THE CALLING PROGRAM 
♦« WITH THE VALUE 1 TO INPtCATK THAT ALL DERIVATIVES SHOULD BE ♦ 

♦ COMPUTED AND Sl’UKEU IN THE FIRST ROW UF THE F ARRAY BEFORE THE^ 

♦ SECUNI) CALL IS MADE UN HAMING.ON SECOND ENTRY TO THE FUNCTIOH ♦ 
♦tDETERMlMED BY LOGICAL VAKIABLE PKEDJ,HAMlNG USES THE HAMING ♦ 

♦ COKKECl’ .K TO COMPUTE MEW SOLUTION ESI IMATES, ESTIMATES THE ♦ 

♦ TRUMCATIJN EKKURS TE(..JJ fur the CURRENT STEP, IMPROVES THE ♦ 

♦CORRECTED SOLUTIONS USING THE NEW rKUNCATlON ERROR ESTIMATES ♦ 
♦SAVES THE IMPROVED SOLUTIONS- IN THE FIRST ROW OF THE Y MATRIX^ 
♦AND RETURNS 10 THE CALLING PROGRAM WITH A VALUE 2 TO INDICATE ♦ 
♦COMPLETl ON • OF OWE INTEGRATinN STEF ♦ 

♦ ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦:♦♦♦♦ ♦♦♦♦♦♦♦ 

INTEGER HAMING 
LOGICAL FRED 

DIMENSION YPReOC2U) ,TE(N) , YC4,N) ,F!:i,N^ 

DATA FRED / , TRUE. / 

C ... 

irC ,NDT,PRSU) GO TO 4 
DU 1 0*1, N ■ 

1 YPREDtOI a Y14,J) + 4 , C 2 , ♦F C I , J J -F t 2 , J ) +2 , ♦F ( 3 , J J ) /3 , 

DU 2 0=1, N - . . . . 

DU 2 R5al,i 
K * 5 " KS 
YU,J)*YCX-lr 0) 

2 irCR,LT.4) FCK,0)«rU-l,J) 

DU ' i 0*1 , N 

3 Y(,l,JJ*YPREUCJ)+112.*rEU^/9. 

XsXfH 

PKEU * .FALSE. 

HAMING * 1 
RETURN 

4 DO S 0*1, N 

YUfJ) * C9,»Yt2,O)-YC4,0)+3,»HVtrCl,0)“l-2,*F(2,0)*FC3,a)))/8, , 

5 Ytl,a)=Y(l,0)-rE(0) 

PEED » .TRUE. 

HAMING * 2 

RETURN 

END 

rUNCTIDiV KUNGKC fJ,Y,F,X,H ) 

C V♦♦V*♦♦♦¥♦’•(♦♦♦♦♦♦♦♦♦¥♦▼♦V♦>)c♦¥♦«♦**4♦V**♦♦♦♦♦♦!V♦♦♦♦♦♦♦ ♦♦♦♦*♦♦♦♦ 

C ♦THE FUNCTliJ!“ KUNGE' EMPLOYS FOURTH ORDER RUNGE-KUTIA METHOD ♦ 

C ♦WITH RUN'JE'S COEFF'S TO INTEGRATE A SYSTEM OF W SIMULTANEOUS ♦ 

C ♦ FIRST ORDER URDIMARX UIFKERENTIAL EGUATIOWS F ( 0 ) *01 C 0 ) /DX ♦ 

C ♦,0*1,N,ACHUSL A5TEP OF LENGTH H IN THE INDEPENDENT VARIABLE H ♦ 

C ♦!« THE INUEPENDEW? VARIABLE X, SUBJECT TO INITIAL CONDITON YCJVt 

C ♦! ,0*1, n, EACH FIJI, THE DERIVATIVE OF YtO),MUST BE CUMPUTED FOOR» 

C ♦TIMES PER integration STEP BY CALLING "THE PROGRAM, THE rUNCT10S» 

C ♦MUST BE CALLEU-FIVE TIMES PER STEP .'SO THAT THE INDEPENDENT ♦ 

C ♦VARIABLE-VALUE (XJ AND THE SULN VALUES ( Y C 1 , I ) . , , Y t U N J ) CAN BE* 

C ‘ ♦OPUATEU USING KUWGE XUTTA ALGORITHM, M-PASS COUNTER,* ♦ 

C ♦♦♦♦*♦*♦♦♦♦♦?♦♦♦♦♦)!(♦♦♦*♦*♦ ♦♦♦♦)((♦♦♦♦♦*♦♦♦!(!♦♦♦ ♦♦*) 

INTEGER KUNGE 

dimension PHICBU) ,SAVEYtSO) ,y cnj ,FCN1 
DATA M/U/ - 
M*M + 1 

GO TO U,2,i,4,5),M 
I KUMGfc*! 

RETURN 
DU 22 J»i,N 


2 



hAVtxuji = nj) 
Fttitujsri. j j 

A»X + U, b»Ml • ■ 

KUNOK s 1 

KfetOKN 

UtJ JJ a = l»N 

I'M! tJ J xlf'M U J ) ,«»!•' U ) 

KUMUKsi - ■ 

Htri’UKN 

yy 94 jBi.N 

XtJJsSAVKrU) +• 

KUNfiieai 

KBTPMKN 
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